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I. INTRODUCTION 



Interpretation of many modern experiments requires a covariant description of the two- 
body system. This is either due to high precision that calls for an inclusion of all possible 
corrections to a standard (possibly nonrelativistic) approach or due to the high energies 
and momenta involved in the processes investigated. Even the simplest atomic object, 
the hydrogen atom, or the simplest leptonic system, the Positronium, needs a covariant 
description to match the high experimental precision already achieved. In the subatomic 
field the most obvious examples are the properties and structure of the deuteron and to 
some extent the mesons, if the later are treated as a quark antiquark system. In the spirit 
of a local quantum field theory the starting point of a relativistic covariant description 
of bound states of two particles is the Bethe-Salpeter (BS) equation. However, despite 
the obvious simplicity of two-body systems, the procedure of solving the BS equation 
encounters difficulties. These are related to singularities and branch points (cuts) of 
the amplitude along the real axis of the relative energy in Minkowski space. Therefore, 
up to now the BS equation including realistic interaction kernels has been solved either 
in Euclidean space within the ladder approximation (see e.g. Refs. and 
references quoted therein) or utilizing additional approximations of the equation itself {(], 
!3oing so, a fairly good description of experimental data has been achieved 



fcf. Refs. 



16|). Note, however, that for processes involving the deuteron the new 



data on electromagnetic form factors [T7| and on the proton-deuteron scattering reactions 
Q are still changing. In addrtron, there are open lsSu e s describrng the po S itroni um a S 
a bound electron-positron state interacting via electromagnetic forces [20( . 

At high energies a consistent relativistic analysis of processes with two particles is of 
a great importance since covariance and relativistic invariance play a crucial role in the 
calculation of matrix elements. At low energies a relativistic framework can also be rele- 
vant, because of the high precision achieved by modern experiments. This means that a 
relativistic framework might be essential to facilitate understanding of the properties of 
bound system even at non relativistic energies. E.g., investigating the energy spectrum 
of the hydrogen atom as a bound system of an electron and a proton that interact via 
the Coulomb potential requires a refined description of the fine structure of the spectrum 
that can be accomplished by using the relativistic Dirac equation for the electron wave 
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function [21(. Otherwise, in a pure non relativistic framework a sizable amount of correc- 
tions to the potential must be considered. A similar scenario appears in hadronic physics 
when analyzing deuteron reactions mentioned above, were the deuteron is treated within 
the impulse approximation. It is common notion that a non relativistic approach fails to 
describe processes at intermediate and high energies unless additional (relativistic) correc- 
tions are taken into account. At the same time, the use of the BS amplitude already in the 
"impulse approximation" provides such corrections, e.g. the Lorentz boost effects, meson 



exchange currents etc 



221 ] , in a natural way. Hence, analyzing different hadronic and 



electromagnetic processes involving light nuclei 13 



14, 



15. 



16, 



22, 



25] we conclude that 



the relativistic description of the initial and final states of the reacting particles should 
be included ab initio. 

A realistic implementation of this program started off by the pioneering works by 
J. Tjon The basic challenge of the BS approach and that of other relativistic 

formalisms, such as the quasipotential approach or equations of relativistic quantum me- 
chanics within Light Front Dynamics, is the covariant description of the nucleon-nucleon 
interactions. As new refined experimental data on deuteron reactions at intermediate and 
high energies have become available (cf. Refs. 17, Q]), the interest in a theoretical 
treatment of relativistic equations has been renewed: the procedure of solving the BS 
equation numerically has been revisited in Ref . 27j , a reduction of the BS equation to an 
equation of the Light Front form has been proposed in Ref. [9|. Also, a detailed investi- 
gation of the solution and the properties of two spinor particle bound state within Light 
Front Dynamics has been reported in Refs. 28|, l29l |. 

Unfortunately, our understanding of the mathematical properties of bound states 
within a relativistic approach is far from being perfect. In mathematical terms the BS 
equation itself is a quite complicated object, and the technical problem of solving it is still 
a fundamental issue. Consequently there are very few successful examples of solving the 
BS equation for fermions including realistic interactions. For example, in Refs. 2], 4] the 
BS equation for spinor particles was regularly treated by using a two-dimensional Gaus- 
sian mesh. That series of studies revealed a high feasibility of the BS approach to describe 
nucleon-nucleon interactions, in particular, processes involving the deuteron. However, 
it should be mentioned that the algorithms exploiting the two-dimensional meshes are 
rather cumbersome and require large computer resources. In addition, the numerical so- 
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lution is obtained as two-dimensional arrays which are quite awkward in practice when 
computing matrix elements and in attempts to establish reliable parameterizations and 
possible analytical continuations of the solution back to Minkowski space. Therefore, it 
is necessary to provide a method for solving BS equations that would feature a smaller 
degree of arbitrariness. 

In the present paper we suggest an efficient and promising method to solve the BS 
equation for fermions involving interaction kernels of one-boson exchange type supple- 
mented by corresponding form factors. It is based on hyperspherical harmonics used to 
expand partial amplitudes and kernels. We show that this novel technique provides many 
insights into the BS approach. The current study is partially stimulated by the results 
reported in Ref. [8j. We explore the structure of 1 S , and 3 Si — 3 Di bound states for dif- 
ferent couplings and study the details of the convergence of solutions and corresponding 
eigenvalues. In particular, on the basis of the proposed method for solving the BS equa- 
tion it becomes possible to analyze the specifics of the problem related to the stability of 
bound states in the BS approach. Besides that, the hyperspherical expansion provides an 
effective parameterization of the amplitude, which is extremely useful in practical calcu- 
lations of observables and in theoretical investigations of the separability of the BS kernel 



with one-boson exchange interaction 30. 13 1 ] . 



Our paper is organized as follows. Section II contains, in the context of scalar particles, 
an overview on the use of the hyperspherical harmonics which displays the basic features 
of the method. In Section III, our generalization of the method for the spinor-spinor BS 
equation is presented. We introduce a new basis of spin-angular harmonics in the spinor 
space and present the appropriate decomposition for the BS amplitude. The correspond- 
ing analytical expressions for the partial amplitudes are found explicitly in Euclidean 
space. Numerical calculations together with an analysis of the stability of bound states 
are presented in Section IV. In this Section the computational algorithm is introduced 
and the results for scalar, pseudoscalar and vector meson exchange kernels are discussed. 
We give our conclusions in Section V. The most cumbersome formulas are collected in the 
Appendix. 
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II. OVERVIEW OF THE METHOD 



The main idea of solving multidimensional integral equations such as the Bethe- 
Salpeter equation consists in finding an appropriate decomposition of the unknown solu- 
tion over a complete basis in momentum space, carry out several integrations over this ba- 
sis analytically and solve the resulting equations with respect to the unknown coefficients 
of such a decomposition. In this way one reduces the dimension of the relevant integrals 
and essentially simplifies the numerical procedure of solving the equation. Usually (see 



e.g 



the corresponding basis is chosen to be the complete set of the two-dimensional 
spherical harmonics Yi m (6,4>) which allows one to eliminate all the angular dependencies 
from the corresponding equation. Since the BS equation is a four-dimensional integral 
equation, the decomposition over the spherical harmonics results in an equation with only 
two-dimensional integrations that has to be solved numerically. Instead, in the present 
paper we suggest a higher dimensional basis, i.e. the basis of hyperspherical harmonics. 
Apart from the familiar spherical angles 6, <p, it also includes a third variable - the hy- 
perangle x ( see ? below). With this basis one can eliminate up to three variables from 
the four- dimensional integral equation and to reduce the problem to find numerical so- 
lutions of ordinary (one-dimensional) integral equations. To demonstrate the essence of 
the method, let us consider the simplest case of the homogenous Bethe-Salpeter equation 
for two scalar particles with equal masses m interacting via exchange of a scalar particle 
271 ]. The corresponding BS equation for the vertex function Q(p) is 



with mass fx 



G(p) = *g 2 I -^v(p,k)s(h)S(k 2 )g(k), 



where k± t 2 = P/2±k, and k = (ko, k),p = {po,p) are the relative 4-momenta, P = (M, 0) 
is the total 4-momentum of the two particles in their center of mass system and g is the 
corresponding coupling constant. In (pQ) S(ki^) and V(p, k) are the free propagators of 
the constituents and of the exchanged particle, respectively 

g(fcl - 2) = fc?, 2 -^ + fe ' v{p > k) = (p-ky 1 -fx* + i e - (2) 

It is seen from fl2j) that even in this simplest scalar case the BS equation contains singu- 
larities (poles) in Minkowski space. Moreover, it is known that the vertex function G(p) 
itself contains cuts along the real axis of po (ko) making the solving procedure rather 
cumbersome. To rid us of difficulties connected to the treatment of these singularities, 



one usually performs the Wick rotation in the complex plane {po,ip4) and solves the BS 
equation ([I]) in Euclidean space, where it is given by 

Q{P4,P)=9 T^TaTZ UT2 I ..2M.2 I ^2 ^2/^2 I ^12 " ( 3 ) 



., (2tt) 4 (p - A;) 2 + fi 2 (k 2 + m 2 - M 2 /4) 2 + M 2 k\ ' 
To reduce the dimension of the integral we decompose the relevant quantities in over 
the basis of hyperspherical harmonics Z n i m (x, 9, 4>) 



Znlm(X,9,4>) = X n i(x)Ylm{0,< 



XM = \j- (w + f + 1)| Bm'xCa(cosx), (4) 

where Yi m (9,(f>) are the familiar spherical harmonics, and C "t * are the Gegenbauer poly- 
nomials. The hyperangle x °f a 4-vector p = (zp4, p) is defined as 

P4 . |p| 

cosx=— ; smx=— (5) 
p p 

with Euclidean 4- vectors p and having modules p = a/p 2 + p 2 and = v/fef-i- k 2 . The 
hyperspherical harmonics (T4]) satisfy the orthonormalization relation 

d9 sinO I dxsin 2 xZ k i m (x^A) Z k'i'm'(x^,(f)) = hk'Su'5 mm '. (6) 



For the interacting kernel in ([3]) it follows 
1 



2ff2 E ^~+i Vn ^' k ) Z nlm{X P , P , <f>p)Klm(Xk, 6k, <f>k), (7) 



(p — k) 2 + u? 

VK ' r nlm 

Vniprk) = (A J a-) 2 (aTTX 



A± = ^/(p±A:) 2 + /i 2 . 
The vertex function is then given in the form 

g(k 4 ,k) = J2<Pi(k)Z nlm (xk,0 k ,(f)k)- (9) 

nlm 

Changing the integration variables, d 4 k = k sin 2 Xk sinOkdkdxkdOkdcpk, inserting (j7|) and 
([9]) into ([3]) and performing the necessary angular integrations we obtain a system of 
one-dimensional integral equations for the expansion coefficients tpf 

r dkk 3 °° 

<p?(p) = / i&rW, *) E ^H^rW- (io) 

^ 771=1 



The explicit expressions for the coefficients Af m (k) result from the corresponding angular 
integrations. Note that a two particle bound state is characterized by definite values of 
angular momenta I, so that only few (one) values of I contribute in (Q and ffTO]) . 

Equation (TTUT) demonstrates how one can obtain an equivalent system of one dimen- 
sional integral equations by starting from the four- dimensional integral equation ([3]) and 
applying the hyperspherical harmonics decomposition ([9]). Formally the expansion that 
leads to the final equation ( TTUT) contains an infinite number of terms, (n,m = l,oo), 
hence, from a pure mathematical point of view, the problem of finding a solution seems 
not much simplified. However, in practise when starting from an approximate solution of 
([3]) , it turns out that the series ([9]) converges rather fast and it suffice to restrict oneself to 
a finite number of terms, making then the system finite and mathematically meaningful. 
Then the procedure of solving the system ( flQj) numerically becomes rather straightfor- 
ward. Previously we found that the first three to four terms in the decomposition (Q 



assure a sufficiently high accuracy of the solution [30|, [31] . Only in the case of a low 
binding energy M and very light exchanged particles, \i — *■ 0, more terms are needed in 
the series (|9]) for a convergence. In the limit /j = 0a more extended analysis is required. 
Note that the described method can also be applied to the inhomogeneous BS equation. 



III. SPINOR-SPINOR BS EQUATION 

Below we generalize the method introduced in the previous section to the spinor-spinor 
BS equation, i.e., the BS equation for two spinor particles interacting via one-boson- 
exchange potentials. In this case the vertex G(p) is a 4 x 4 matrix in spinor space and 



the corresponding BS equation is [32] 



Sip) = ^ / (01 V <P> k ) r W W S{k) S(k 2 ) f(2), (11) 

where the propagator V(p,k) for the exchange of scalar and pseudoscalar mesons is given 
in (J2J) and by 



■r 

q £ — fi z + is 



V(p,k) = - f-fL- (12) 



for vector mesons. The propagators of the spinor constituents are 

s(h) = 12 kl \ m . , §(k 2 ) ee cs(k 2 fc h - "' 



kf — m 2 + is ' " k\ — m 2 + is ' 
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with the charge conjugation matrix C = ry°7 2 . The meson vertices V are determined 
by the corresponding effective interaction Lagrangians describing the interaction of the 
spinor particles with mesons. For a system of two nucleons they are given by 

T(l) = 1 f (2) = -1 scalar (13) 

T(l) = 75 T(2) = —75 pseudoscalar (14) 

r(l) = 7m - ^ °» P q P f (2) = lv + ^- a up q" vector (15) 

In f fT5|) the momentum transfer q is defined as q = p — k, k denotes the strength of the 
tensor part of the interaction and the coupling constant g in (TTTT) is imaginary for the 
pseudoscalar mesons else real. Each interaction vertex T is regularized by a monopole 
form factor 

A 2 - u 2 

m - ^ (i6) 

with A as free cut-off parameters. For the sake of simplicity we presently adopt k = 0. 
Consequently, in the propagators of vector particles, eq. ( [121) . the part proportional to 
q^qv/fj* 2 is neglected as well. These restrictions lead only to slight redefinitions of other 
effective constants (cut-off parameters, vector coupling constants etc) and do not affect 
the method and the main final conclusions. 

Contrary to the scalar case, the bound state of two spinor particles is characterized by 
the total angular momentum J, which is an algebraic sum of the total spin S = 0, 1 and 
total orbital momentum L (L = 0, 1,2....), i.e., J = S + L. Traditionally, for two spinor 
bound states one adopts the spectroscopic notation 2S+1 Lj. 

We investigate the lowest ground states of the 1 Sq and 3 S\- 3 Di channels taking into 
account only one type of mesons at a time, either scalar, pseudoscalar or vector. The 
inclusion of the sum of all mesons in the interaction kernels of the BS equation which 
reflects the more realistic case is beyond the scope of the present paper. A generalization 
of the method will be done elsewhere. 

For completeness we give the BS amplitude ^ used in the following that is related to 
the corresponding vertex functions via 

*(p) = s( Pl )g( P )s(p 2 ). (17) 



8 



A. Spin-angular harmonics 



The main difference between the scalar ([3]) and the spinor BS equation (TIT!) is that in 
the later case the vertex function Q(p) is a 4 x 4 matrix in spinor space. Consequently, the 
spinor BS equation is of matrix form and before proceeding with numerical solutions we 
shall transform it into a system of ordinary equations. To this end, we expand G(p) into 
a complete set of 4 x 4 basis matrices and obtain a system of equations for the coefficients 
of such a decomposition. In the most general case, there are 16 linearly independent 



matrices that can be used as a basis. The choice of such a basis system depends on t 
specific aim envisaged. Usually one uses either the complete set of the Dirac matrices 



he 

41 

26j. Different basis are related to each 



or the complete set of the spin-angular harmonics 
other via unitary transformations [ljj. 

For specific bound states with given quantum numbers only some basis matrices con- 
tribute to the vertex function Q(p) (amplitude ^(p) ). E.g., for the 1 Sq state only four 
matrices are relevant to describe the amplitude, while in the 3 Si — 3 D\ channel eight 
basis matrices are needed. A standard choice for the basis matrices in these two channels 
is referred to as p-spin angular harmonics T a (p), where the index a lists the quantum 
numbers of the LSJ momenta and that of the p-spin. In fact, the spin angular harmon- 
ics are constructed as outer products of two Dirac spinors which, for each constituent, 
form complete sets of solutions of the free Dirac equation. Then the p-spin (projection) 
reflects the signs of the energy of two basis spinors in the outer product and is labeled by 

++, ,H — and — h correspondingly. Sometimes, instead of H — and — h one uses linear 

combinations which can be even (e) or odd (o) with respect to symmetry on the relative 
energy po (for details see Ref. [26]). In Table [T] we present the classification of the par- 
tial components in this basis and the spectroscopical notation for the partial components 
of the BS amplitude in the 1 Sq channel (first row) and in the 3 Si- 3 Di channel (second 
row). The explicit expressions for the corresponding spin angular harmonics T a (p) can 
be found, e.g., in Ref. 

The expansion of the BS amplitude *&(po,p) and of the vertex function G(po,p) into 
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spin-angular harmonics is 

*(po,p) = J> a (po>|)r a ( P ), (is) 

a 

G(P0,P) = ^29a(P0,\p\)Ta(-p), (19) 

a 

where the expansion coefficients ipa{po, \p\) and g a {po, \p\) are evaluated numerically. A 
comparison of (fT9j) with (jlj) and ([9]) shows that the spin-angular harmonics T a can be 
considered as a direct generalization of the spherical harmonics Y\ m in spinor space. Hence, 
at first glance, it seems natural to expand also the coefficients g a (Po, \p\) (^(Po, \p\)) into 
Gegenbauer polynomials X n i(x) defined in (jlj) and to obtain a system of one dimensional 
integral equations like ( ITUl) . However, a more detailed inspection of the spinor-spinor 
BS equation, together with the explicit forms of the spin-angular harmonics shows that 
the use of the commonly accepted form of the p-spin basis hinders a further use of the 
Gegenbauer polynomials. To be more specific, note that in (fT8j) and fjT9l) terms appear that 
always mix the angular and hyperangular dependence (viz. nonlinear terms proportional 
to, e.g. l/|p|) which make an analytical integration over the hyperangle x impossible. To 
avoid such difficulties within the hyperspherical harmonics formalism we suggest a slightly 
modified set of spin-angular matrices which represent a generalization of the p-spin basis. 
In the 1 5'o channel the new basis is 





l 

Vievr 75 ' 


Up) = 


vT6^ 7075 ' 


Up) = 


1 (p7) 


V167T |p| 


Up) = 


1 (pj) 


V16tt |p| 



7o75, 
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and for the 3 S\- 3 Di channel we choose 



Tib) 

Up) 
Up) 
Up) 
Up) 
Up) 
Up) 
Up) 



i 



1 



(7* 



M)i 



7o(7^m)> 



VT6^ 

VTfo \p\ ' 
^3 1 7o 



2^/16% \p 
V2 1 



2 v 7 !^ 
\/2 1 



[(P£ M ) + (P7) (7^)] > 
3 



IP 



(211 



2 Vi6 



=7o 



7T 



|P| 



(p£m)(P7) 



7o- 



^3 1 



|P| 
1 



[(p£m) + (P7) (7^)] • 



2 vT6^ IpI 

The left hand side of (121!) depends implicitly on .M, which denote the components of the 
polarization vector £ M fixed by £ +1 = — (1, i, 0)/\/2, = (1, — 0)/\/2, £ = (0,0,1). 
The new basis is orthogonal and normalized 



(22) 



(23) 



dVt p Tr [^(^^(p)] = 5 TOn 5xM'- 
The partial expansions of the vertex functions over the new basis is given by 

G(po,P) = ^29n(po,\p\) T n(-p)i 

n 

with 

9n(Po, |p|) = J dQpTr [g(p ,p)T+(-p)]. 

In (1221) n — 1 . . . n max , where n max = 4 for the 1 5o and n max = 8 for the 3 5i — 3 D\ 
bound states. As mentioned, different bases are related to each other via corresponding 
unitary transformations. The connection between the p-spin and the new basis (|21|) . or 
equivalently, the relation between the quantities g a (a = JLSp) and g n (n = 1 . . . n max ) 
can be found by using the completeness of the two basis and the parity of the coefficients 
Qi with respect to the operation po — > —po- The Pauli principle together with charge 
conjugation operation C leads to 



cr n T (-p)c 



7T 



1 



J-l 



Up), 



(24) 



n 



where I stands for the isospin of the system, tt Po is the p-spin parity. From (I24I) one 
obtains that in the 1 S'o channel the component g^ is of the odd parity (tt Po = — 1) while 
gi, ...,#3 are of the even parity (tt po = +1); correspondingly, in the 3 Si — 3 D x channel the 
two components #7, #8 are °dd in p , and gi, ...,ge are even. 

The unitary relation between two sets eq. (TT9~j) and eqs. (12"U|) and (12"T1) is 

IpI 



9^8++ 



1 1 m 

7! 91 + Vl^ 92 " 

1 1 m 



m 

ti p ti p 



gap? 



-93, 

-93, 

) 

94, 



(25) 



for the 1 Sc) states and 



9^8++ 






/ y/2m+2Ep 
6 E p 


2m+E p 
6 E p 


i_m 

s/E Ep 


1 IPI 

Ep 


1 Ep—m 
3 Ep 


1 Ep—m 
3 Ep 


\ 


M 








y/2 m+2E p 


y/2 2m+E p 


i_m 


1 IPI 


1 Ep—m 


1 Ep—m 




^2 






6 E p 


6 E p 


V6 Ep 


Ep 


3 Ep 


3 Ep 




9s D ++ 




1 Ep—m 
3 Ep 


1 Ep—m 
3 Ep 


1 IPI 

V3 Ep 


i_m 

s/E Ep 


y/2 2m+E v 
6 E p 


y/2 m+2Ep 
6 Ep 




^3 


9*s— 






1 Ep—m 
3 Ep 


1 Ep—m 
3 Ep 


1 IPI 

Ep 


J_]PI 


2m+Ep 
6 Ep 


m+2Ep 
6 E p 




94 


9^P{ 






1 IPI 

V^3 E p 





I m 
Ep 





3 Ep 







95 


\ 93 pi 


) 




V 


3 E p 





m 
Ep 





1 IPI 

%/3 Ep 


/ 


w 



(26) 



for the 3 S*i — 3 state. Here E p denotes the total energy of a free particle with the 
momentum p, i.e. E p = p 2 + m 2 . As seen from (I25I) and (I26I) all quantities containing 
Ep and |p| that prevent the use of Gegenbauer polynomials have been explicitly extracted 
into the corresponding coefficients. 

To complete this section we present the partial amplitudes in terms of the corresponding 
partial components of the vertex functions in the p-spin basis 

9s(D)++(Po, |p|) 



^S(D)++(P0, |P|) 
$S(D)—(P0, |P|) 

i>Pe{po, IpI) 
i^Poipo, IpI) 



(M/2 

9s(d)—(pq, IpI) 

(M/2 + E p 
g Pe (Po,\p\)(M 2 /4 



2 ' 
Po 



2 ' 
Po 



Pi ~ E p) + 2 9p (Po, \p\)PoE p 



(M 2 /A-p 2 -E 2 ) 2 -Ap 2 E 2 
9p (Po, |p|)(M 2 /4 -vl ~ E 2 p ) + 2g Pe (p , \p\)p E p 



(M 2 /A - p 2 



Elf 



(27) 
(28) 
(29) 
(30) 
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Since the components with negative p-spins reflect the contribution of negative energies 
(of the solution of the Dirac equation) one might argue that at moderate energies and 
momentum transfers they could be suppressed compared to the components with positive 
p-spins. Hence in calculations of matrix elements such components could therefore safely 
be neglected. Moreover, often the components with mixed p-spins (+— or — h) could also 
be disregarded in concrete calculations. Note, however, that in the BS equation these 
components cannot be omitted, since the convergence and stability of the numerical solu- 
tion is sensitive to each component. To estimate the magnitude of different components 
one might introduce the pseudo probability, i.e. the contributions of each component in 
the normalization condition 

P S (D) ++ = J dpod\p\(E p -M/2)\^ s{D)++ (p ,\ P \)\ 2 , 
Ps(D)- = - j dp d\p\(E p + M/2)|V5(d)-(Po>|)| 2 , (31) 
Pe = -M J dp d\p\\i;p e (p ,\p\)\ 2 , 



Po = -M J dp d\p\\ip Po (p ,\p\ 



B. Hyperspherical decomposition 



Similar to the scalar case the spinor-spinor BS equation after the Wick rotation [36], 
is considered in Euclidean space (c.f. eq. Q). The Wick rotation can be achieved 
by replacing the scalar product of two vectors in Minkowski space by their Euclidean 
analogue (k 2 = h% — k 2 — > k 2 = k\ + k 2 ) and changing p — > ip 4 and k — > ik 4 . Notice, 
that every odd function of po, being homogeneous in its argument, explicitly receives an 
additional imaginary unit i after the Wick rotation. Therefore, for convenience, we divide 
this common factor by redefining the odd partial components 

for the 1 S'o channel and 

97 -> 197, 9s -> m- 
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for the 3 Si — 3 D\ channel. Now, by placing eq. fl22|) into the Wick rotated BS equation 
( fTTj) and using fl23l) one obtains 

9n(P4, \p\) = 9 2 b n J dVL p J ^^ I 5 , (/c4, |fc|) ^ - fc j 2 - ^ 2 ^2 A nm(p, k)g m (k A) k), (32) 

where m, n = 1 . . . 4 for the 1 S'o and m,n = 1 ... 8 for the 3 S\- 3 Di state, respectively, b n 
is a sign coefficient reflecting the type of the exchanged meson (see Tables [III and IHIj) and 
the scalar part S(k^, \ k\) of the two spinor propagators is defined as 

S(h,\k\) = — 5 • (33) 

(P + m2-Mi) 2 + M2fc 4 2 

The coefficients A nm (p, k) in (1321 directly follow from calculating traces of the BS equation 
after multiplying it from the right by the corresponding basis matrix (120]) or (T2~T|) . The 
angular dependence of A nm (p,k) on f2 p is entire ly g overned by the dependence on the 
vector p of the corresponding harmonics T n (p) [421 ] . So, from ( l20i) one infers that in 
the 1 S channel A Xm and A 2m are oc Y 00 (p) while A 3m and A im are oc ^-^(fc)^!^^). 
Analogously from eq. (I2T!) it can be found that in the 3 Si — 3 channel the angular 
dependence of A nm (p, k) is proportional to Yi^{p) with / = 0,1,2. Such a dependence 
essentially simplifies integrations over Q p . 

After expanding the interaction kernel into hyperspherical harmonics, eq. ([7]), integra- 
tions over the angles Q p and Qk are carried out analytically. The result is 

9n(P4,\P\) = 9 b n t-t^tt S (k A , \k\)W ln (p, k , (h, k)g m (k 4 , k), 

(34) 



where 



Wi n (p, k, x P , Xk) = Jl~i Vl ^ k)Xu n (x P )Xu n (xk 



i + 

is the hyperspherical partial kernel in Euclidean space. The values of the angular mo- 
mentum l n are restricted by the dependence of A nm on Q p of the corresponding integrals. 
These are l n = for n = 1, 2, and l n = 1 for n = 3, 4 in the 1 S'o-channel and l n = for 
n = 1, 2, /„ = 1 for n = 3, 4, 7, 8 and l n = 2 for n = 5, 6 in the 3 Si — 3 -Di channel (see 
also eqs. ( 1201 and ( 12TI) ). The explicit expressions for the quantities a nm in eq. ( 1341) are 
collected in Tables HVl IVl and I VII A prominent feature of the partial BS equation (l34"l) that 
is related to the proposed basis ( 120]) and ( |2T|) is that the coefficients a nm are independent 
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on the type of the exchanged meson in the interaction kernel. This dependence is solely 
contained in the coefficients b n (Tables [Til and III II ) and in the coupling constant g (if a 
cut off form factor is included into the consideration, the cut off parameter A can also 
depend on the meson type). From Tables UTtiVIl one can also infer the parity tt po of the 
component g n (j)4, \p\). Note that in eq. (IMI) one can easily express S , (fc 4 , |fe|) of eq. (1551) 
and a nm (k^, \ k\), given in Tables HV ] IVl and IVI | via the hyperspherical variables Xk and k. 

We now expand the partial vertex functions into hyperspherical functions in a similar 
way done in eqs. ([7]) and (jHJ). Since the value of the angular momentum l n are restricted 
in the partial components g n , as discussed before, the hyperspherical expansion reduces to 
an expansion into the functions Xji(x p ), i.e. the Gegenbauer polynomials. Moreover, due 
to definite parity of the components g n with respect to cosx P (cosx P = Pa/p), summation 
over j is restricted to only even or only odd values of j, in accordance with the relation 

Cj(-cos Xp ) = (-iy'Cj(cos x P ). 

For the partial BS components g n one gets 

oo 

01,2(P4,|P|) = $^^l 2 (p)^2i-2,o(Xp), (35) 

oo 

93(Pa,\p\) = ^2gi(P) X 2j-i,i(x P ), (36) 

3=1 
oo 

9a(pa, |p|) = ^2gi(p)X 2j ,i(x P )- (37) 
i=i 

Similarly for the 3 S\- 3 Di channel we obtain 

oo 

9ij(P4,\p\) = ^29i 2 (p) X 2j-2,o(x P ), (38) 

3=1 
oo 

93,i(pi,\p\) = ^gi,i(p)X2j-i,x{Xp), (39) 

oo 
oo 

97,S(P*AP\) = 5Z#7,8(p) X 2i,l(X P ), ( 41 ) 

where the corresponding equations for the coefficients g^ip) can be readily obtained by 
inserting eqs. (|35l) - (|41|) into eq. (|34|) and by using the orthonormalization relation for 
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Xji(Xp)> e 1- &■ m the 1 Sq channel the resulting equations are 



01,2 (p) 

5a (p) 
04 (p) 



-g 2 b 



1.2 



dA; fc 3 



4 oo 



^ 2 (2j 



n=l m=l 



"A 



dfc A; 3 
8tt 2 2j 



4 oo 



dk k 3 



n=l m=l 

4 oo 



8tt 2 (2j + 



(42) 
(43) 
(44) 



n=l m=l 



For the 3 Si- 3 Di channel we get 



01,2 (£) 

03,4 (?) 
05-8 (p) 



-0 2 & 



1.2 



8 oo 



?7T 2 (2j - 1) 



n=l m=l 



-9 2 b: 



'3,1 



dA; A; 3 



8vr 2 2j 



dfc A; 3 



n=l m=l 

8 oo 



8vr 2 (2j + 1) 



(45) 
(46) 
(47) 



n=l m=l 



where the explicit expressions for A^(fc) and BJ^(h) are presented in Tables IVHl IVIIII 
and [IX] (for details see Appendix |A]). This is our main analytical result. As seen from 
Tables HT1 and lITlj the partial amplitudes g% in the 1 Sq channel and g2, g% and g$ in the 
3 Si — 3 D\ channel are identically zero for the vector meson exchange, which is a direct 
consequence of our (approximate) choice of the propagator of vector particles. 

An implementation of the vertex form factors (fl~6l) into the calculation is straightfor- 
ward. Namely, by observing that in Euclidean space 

(p _ fc | 2 + /i2 F[(p - k) 2 ] 2 = 2tt 2 ]T -1- v n (p, k) Z nlm (x P , P , 4> P ) Z* lm ( X k, k , <f>k), (48) 

with 



nlm 



Vn(P,~k) 



(A£ - At) n (Ai - A A ) n 



(A^ + A^)™+ 2 (A^ + A A )™+ 2 



a( i n A2 ~^ 2 ( A +~ A -) n 
{n+ ' A A A A (A£ + A*)"+ 2! 



(49) 



A^ = a/(p±*0 2 + pA A^ = V(P±^) 2 + A 2 
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one finds that the only modification of eqs. (142]) - (14?]) consists in replacing the quantities 
V n (p, k) of ([8]), by the "dressed" quantities V n (p, k) of (|49l) . Hence, one can conclude that 
the system of equations (14"2"1) - (1471) together with the proposed new basis (T2U1) and (T2"TT) 
present a direct generalization of the hyperspherical harmonics method for the spinor- 
spinor BS equation. Since the new basis used here is less intuitive than the more tradi- 
tional p-spin basis, one would like to relate the results to the p-spin formalism where the 
computed matrix elements have a clear physical meaning and allow for direct comparisons 
with nonrelativistic calculations. This can be done by first solving the BS equation for the 
partial components g n and then by use of eqs. ( 1251) and ( 1261) finding the desired partial 
components g a of the p-spin formalism. 



IV. NUMERICAL SOLUTIONS 



Equations (H2l)-( [4T1) represent the desired system of the BS one- dimensional integral 
equations within the hyperspherical harmonics formalism to be solved numerically. Before 
choosing a specific computational algorithm one has to analyze at least two issues, i) 
existence and uniqueness of the solution and ii) if such a solution exists, one has to 
analyze the convergence and stability of the (approximate) solution to the exact one. 
Obviously, both issue are tightly connected to the properties of the interaction kernels. 
The main requirement for the existence of solutions of the Fredholm type equations is the 
finiteness of the kernel. This can always be fulfilled by considering the cut-off form factors 
(1T6|) in the interaction Lagrangians. To find the numerical solution, one transforms the 
continuous space of arguments and solutions into a discrete one, forms the skeletons of the 
approximate kernel and of the solution, chooses a numerical method to solve the resulting 
(finite) system of linear equations and investigates, within such a scheme, the convergence 
of the skeletons to their exact originals. In such a way one estimates the effectiveness and 
correctness of the algorithm of the chosen procedure. However, inclusion of cut-off form 
factors is not always physically justified. Moreover, in some specific cases the interaction 
kernel can increase unbound and therefore, does not automatically guarantee the existence 
of a solution. Within nonrelativistic quantum mechanics it is known that the bound state 
of the Schrodinger equation with the interaction potential of the form U(r) = —a/r 2 
(r is the radius vector in co-ordinate space) can disappear [43| at large enough coupling 
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constant a 



2l| . A similar situation occurs also in relativistic quantum mechanics within 



the Light Front Dynamics 



37| |. where for potentials corresponding to exchanges of 



scalar mesons, a critical value g CT of the coupling constant exists, above which the bound 
state disappears. Recall that this refers to the case without cut-off form factors. Inclusion 
of form factors, e.g. eq. (1161) . essentially aims to improve convergence of the system of 
equations and to assure the existence of a solution for any type of the (attractive) exchange 
meson. 



A. Stability of solutions 

In this section we investigate the existence and stability of the solution of the BS 
equation without cut-off form factors for two spinor particles. For defmiteness, we consider 
the system ( 132]) in the 1 So channel for the scalar exchange. 

Since we are interested in clarifying the existence of a critical value of the coupling 
constant g CT , the system is analyzed at asymptotically large values of p. In this 
case it is sufficient to investigate the properties of the system (1321) in the two-dimensional 
space (jp,Xp)i without further decomposition (135]) - (135]) into Gegenbauer polynomials. The 
equation for the component g\ is then 

/- \ __ g 2 f d ~ k k 3 dxk sin 2 \k dVt k 1 ^ 

9i(P,X P ) (2n) 4 J (k 2 + m 2 - M 2 / A) 2 + M 2 k 2 cos Xk (P~ k ) 2 + ^ 2 
x [(M 2 /4 + m 2 + k 2 ) gi {k, X k) + Mmg 2 (k, X k) - Mk sin X k g 3 (k, X k)]- 

For further convenience we introduce a new integration variable 7 by k = 7p. Then at 
p — > 00 in eq. ( 150]) only terms proportional to k survive and eq. ( 150]) becomes 

OO 7T 1 

gi(p,Xp) = J dl j sin2xkdXk J dcos °kV(j, cos 9 k ,x P , Xk)gi(iP,Xk), (51) 

-1 

where the interaction kernel 

V(7, cos fc , x P , Xk) = -7- — 2 — o ~ — o T~- : ' ( 52 ) 

1 + 7^-27 cos x P cos Xk - 27 cos k sin x P sin Xk 

is positively defined and the relation ^(7, cos^fc, Xp, Xk) = cos 8k, Xp,Xk) holds. 

From eqs. (1511) and (1521) one infers that the integral converges if g\ (p, x P ) asymptotically 

vanishes as 1/p or faster, i.e. if g\ can be written in the form 

9i[P,Xp) = ~ 1+0 ■ , (53) 
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where sin Xp has been introduced for convenience. Now, by splitting the integration over 
7 in flSTl) into two sub-ranges as 



d 1= / . . . c?7 + ...d~f, (54) 



and changing the variable in the second integral 7 — > I/7 and carrying out integration 
over cos 9k analytically, we obtain for (l5Tj) 



1 7T 

/— / cosh (/5 In 7) In ^ QCOS (x P + Xfc) fr(fc, Xfc ) = jjZL^^p x ) (55) 

J J J 1 - acos(x P - Xfc) r 







where a = - — ; — -. Eq. (Tooj) can be considered as an equation of the Sturm-Liouville-like 



_2t_ 

I + 7 2 

(27r) 3 

problem of finding the eigenvalues — — and the eigenvectors h(p, x P ) of the corresponding 

9 2 

integral operator. It is immediately seen that the eigenvalues depend on the asymptotic 
behavior of g\ and that the most harmful situation occurs at (3 = and h(p, x p ) = h(x P ). 
Namely, for j3 < the integral is divergent and the equation becomes meaningless, for 
(3 > the integral converges rather fast ensuring the existence of bound states. Thereby, 
the critical value of the coupling constant g cr can occur only at /3 — 0. Note, that from 
eq. (155]) one finds that the function h(p, x P ) is odd with respect to the variable x P i 
h(p,—Xp) = —h(p,Xp)- This implies that at f3 = and h(p,Xp) = h(x P ) the quantity 
h(x p ) can be developed into an odd Fourier series 



KXp) = ^c n sm(nx P )- (56) 

n=l 

It is straightforward to check that each term of the series (156]) represents a solution of the 
eq. (155|) with /3 = 0. Moreover, since 



1 7T 1 2vr 

— / «Xfc In- -. r sm(nxfc) = / — / df-- cos(n(£ - x P )) 

7 J 1 - acos(x P - Xk) J 1 J nl-acos4 

00 00 

= i / *I R ( e -*» ) = * S i„(^ p ), (57) 

n J 7 y 1 — a cos 4 rr 

I J 

the critical value <^r depends on n, being g cr = 27rn. Obviously, the lowest g cr occurs at 
n = 1, i.e. # cr = 2tt. 
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Similarly one finds the critical value of the coupling constant from the asymptotic 
equation for g?,. The result is that the lowest critical value is determined by the solution 
of the form h(x p ) = sinx P + esin3x p (with e = 1.100925), which provides g 2 T ~ 71.71862. 
The remaining components and g± are negligibly small and can be neglected in the 
present analysis. Hence, we find that in the 1 Sq channel the critical value g 2 r ~ 4tt 2 . This is 
confirmed, with good accuracy, by concrete numerical calculations where no bound states 
for g > g CT are found at p —>■ oo. This is also in full accordance with the nonrelativistic 
potential V(r) ~ —a/r 2 . 

The critical value for the coupling constant for the 3 Si — 3 D\ channel above which the 
bound state disappears is g 2 r ~ 78. 



B. Numerical methods 



In concrete numerical calculations we form the skeletons of approximate solutions and 
kernels by using the Gaussian method of computing integrals and by restricting the in- 
finite sum over m in eqs. (l42l) - ([44l) and (H5l) - (H71) by a finite value M max . The Gaussian 
quadrature formula assures a rather good convergence of the numerical procedure and 
provides the sought solution in the Gaussian nodes which are spread rather uniformly 
in the interval < p < oo. In order to have the solution in detail at moderate values 
of p, which is the interval of the actual physical interests, one usually redistributes the 
Gaussian mesh making the nodes more dense at low values of p. To this end one applies 



an appropriate mapping of the Gaussian mesh by changing of variables as [27 ] 



1 + x 

p = p(x) = c (58) 

1 — x 

with Cq ctS cl free parameter and — 1 < x < 1. Then the corresponding set of linear 
equations reads as 

X = g 2 AX, (59) 

where the vector 

X T = ( [W(*0}Sfe • • • > [{9n(ki)}^r) (60) 



represents the sought solution in the form of a group of sets of partial wave components 
9ni n — 1, ■ ■ ■ , n max ; m = 1, . . . , M max specified on the integration mesh of the order N G . 
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As before, in eq. (160]) n max = 4 and n max = 8 for the 1 S'o and 3 Si — 3 -Di channels, respec- 
tively. The matrix A is determined by the corresponding partial kernels, the Gaussian 
weights and the Jacobian of the transformation (I58p and is of the N x N dimension, 
where N = n max x M max x N G . Since the system of equations (J59|) is homogeneous, the 
eigenvalues g 2 (at given mass of the bound state M) are obtained from the condition 
det(g 2 A — 1) = 0. Then the partial components g™ are found by solving numerically the 
system ( 159]) with these eigenvalues g 2 . 



We use a combined method of finding the solution X. First, the Gauss- 



ination and pivoting method involving the choice of the leading element 38 



ordan elim- 
is applied. 



Then the obtained solution is used as a trial input into an iteration procedure to find 
(after 5 — 10 iterations) more refined results. Within such an algorithm we investigated 
the convergence of the approximate solution by increasing the dimension of the matrix 
A and found that the method is stable and robust up to values N ~ 4000. This is quite 
enough to obtain solutions with practically any desired accuracy. 

Another way to solve the system (159]) consists in directly using the iteration method 
with known nonrelativistic wave functions as trial inputs. This method has been inspired 
by the success of the One Iteration Approximation scheme developed in 14j, [l5j, [l6j , which 
provides a quick and accurate solution. In the present paper both of these methods are 
widely explored. In practice, however, instead of finding the eigenvalues g 2 (M) at given 
M, one usually considers the inverse problem when the coupling constant g 2 is kept fixed 
and the mass M(g 2 ) is assumed as a function of g 2 . 

As an illustration of the stability of the numerical procedure, in Table |X] we present 
results for the masses of the bound state M(g 2 ) depending on the Gaussian mesh N G 
and M max . Calculations have been performed for the 1 S'o state, eqs. (S2J)-(j44j), with a 
scalar meson exchange of mass fx for two values fx = 0.15 GeV/c 2 and fx = 0.5 GeV/c 2 ; 
the constituent particles (nucleons) have been taken with equal masses m = 1.0 GeV/c 2 
for simplicity. Results presented in Table |X] clearly demonstrate that the approximate 
solution converges rather rapidly, and already at M max ~ 4 — 5 and N G = 64 the method 
provides a good solution of the system (|42p - (j44p . Obviously, the coupling constant g 2 
must be taken small enough, i.e. g 2 < g 2 r ~ 40 to ensure the existence of the solution. 
The free parameter cq in (158]) does not affect the convergence and for definiteness it has 
been set Co = 1. As expected, an implementation of the cut off form factors (TT61 . ( 14*81) 
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and (14*9]) essentially improves the convergence of the approximate solution. 

To display the behavior of the components of the vertex function, in Fig. [1] we 
present the resulting coefficients for the 1 S channel g{(p), j = 1...4, eq. (I3"5"j) . at 
M = 1.937 GeV/c 2 for N G = 96, M max = 4 and g 2 = 15. It can be seen that at large 
p these functions decrease as inverse powers of p which permits to cast the approximate 
solution into the form 



g{{p) 



p 2 V' 1 ^ (61) 



[P 2 + b 2 \ j^(P 2 + b 2 ) 

where the parameters bj and a l j can be found from an x 2 analysis of the adjustment of 
eq. (1611) to the approximate solution (see Table IXlT) . The solid lines in Fig. [T] reflect the 
result of the fit of the numerical solution by eq. ([6"Tj) . The accuracy of the results implies 
that in such a way one can find solutions of the corresponding system of equations as 
continuous functions of p which are extremely useful in practical applications. Similar 
analysis of other coefficient functions g^s^ip) shows that an excellent fit of the numerical 
solutions can be achieved for all the partial components g a {P)- From this encouraging 
result we argue that this method using hyperspherical harmonics can be considered as a 
reliable tool to solve the Bethe-Salpeter equation numerically, even with different param- 
eterizations of the solution than the simple form (IBTI) . Unfortunately, for small meson 
masses \i ~ the method becomes less effective, having a poor or even failing conver gen ce 
feature and hence becomes less adequate, requiring a separate analysis (see also Ref. |27|). 

Also, as mentioned above, at large values of the coupling constants g 2 , close to or even 
larger than the critical value g 2 r , the numerical solutions become strongly dependent on 
M m3jX , Nq and Cq which is a clear signal that (in absence of the cut-off form factors ([161) ) 
the solution becomes unstable at g 2 ~ g 2 T and disappears at g 2 > g 2 T . Such a situation is 
illustrated in Fig. [2] where we present the dependence of the mass M on the values of the 
cut-off parameter A at different coupling constants g 2 , below and above the critical value 
g 2 r ~ 40. It is seen that if the solution exists and is stable (g 2 < g 2 r ) then the mass M 
is almost independent on the cut-off form factor and tends to a constant value at large 
A. Contrarily, the dependence of the mass M on A at g 2 > g 2 r evidently indicates that in 
this case the solution becomes unstable and at large values of A, A — > oo, it can disappear 
at all. Such a behavior of the solution at g 2 ~ g 2 r exactly reproduces the peculiarities of 
the well-known collapse phenomenon for potentials like —a/r 2 in nonrelativistic quantum 
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mechanics (see, also discussions in Ref. 



281]). In Fig. [3] the numerical solutions for the 



partial components g{,j = 1, 2, 3 are given for two values of the coupling constant, below 
(solid lines) and above the critical g CT (dashed lines). In order to ensure the existence 
of the solution, calculations have been performed at a finite, however large, value of the 
cut-off parameter A = 500 GeV/c. It is seen that in the case of g 2 > g 2 r the asymptotic 
decrease of the wave functions is rather weak which implies the instability of the solution. 
In Fig. [3] the continuous lines reflect the result of a fit by eq. fjoTl) which for the asymptotic 
region of the component g\, can be written in a simple form 

MP) - (62) 

where c = 0.77 (g < g cr ) and c = 0.4 (g > g CT ). A comparison with eqs. (1531 and (1561) 
shows that for coupling constants below the critical value, g < g cr , one has (3 > 0, while 
above g < g CI , (3 < 0, i.e., the integral (15T1) without form factors (fl6l) diverges. A similar 
behavior of the numerical solution occurs in the 3 Si- 3 D\ channel as well (see, Fig. HI). 



C. Scalar coupling 

Having established the main features of the numerical procedure, we solve the BS 
equation within the hyperspherical harmonics method for different kinds of the exchanged 
mesons, i.e. scalar, pseudoscalar and vector mesons. We investigate the dependence of 
the bound state mass M on the coupling constant g and study the partial components 
ga as a function of the hypervariable p at fixed values of M and the coupling constant 
g. All calculations have been done with and without the cut-off form factors (flB]) and, 
obviously, for the coupling constants below their critical values. We compare our results 
with other calculations performed for similar conditions, namely we widely compare our 
analysis to the ones obtained by Light Front (LF) dynamics Q] and by the non relativistic 
(NR) Schrodinger approach. Note that within these approaches the dynamical variables 
used differ slightly from the ones used in BS formalism, hence a direct comparison of the 
calculated quantities is hampered. However, one can reconcile approaches by choosing 
one variable within the BS formalism, e.g. the modulus of the 3-dimensional relative 
momentum p, and relate the corresponding quantities in LF or NR approaches through 
this variable for further comparisons. 
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Such relations have been found and reported in detail in Refs. [24], |40(. Here we note 
only that in determining relations between the BS and LF amplitudes one finds that in 
the 1 5'o channel only two components ( 1 5' ( ^ + and 3 -Pq) of the BS amplitude correspond 
pairwise to the two LF functions f\ and f'2 (notation as in Ref. fl). The other BS com- 
ponents, when being projected on to the LF surface, vanish. Similarly, in the 3 Si — 3 D\ 
channel only five BS components survive on to LF surface {24]. ifcl]. In Fig. [5] the mass of 
the bound state as a function of the coupling constant for scalar meson exchange is shown. 
The solid line corresponds to results within the BS approach, while the dashed and dotted 



lines show calculations with LF |8| and NR (with an Yukawa-type potential) approaches, 
respectively. Only at low values of the binding energy different approaches provide similar 
results. As the coupling constant increases (increase of the binding energy) the difference 
becomes more and more significant. Even within two relativistic frameworks, BS and LF 
formalisms, the difference increases with the binding energy increase. This is illustrated 
also in Figs. [6] and [7] where the wave function eq. ( |27l) . is compared to the corre- 

sponding LF wave function f\ [8j for two different binding energies. It is seen that the 
two approaches provide similar results for low, and rather different results for high values 
of the binding energy. Obviously, this is a direct consequence of the different treatment 
of the relativistic effects within the BS and LF formalisms. Since when projecting the BS 
amplitude on to the LF surface, some partial components (with negative p-spins) vanish, 
the observed large difference between BS and LF results serves as a hint that in this 
case the role of the components with negative p-spins increases. This can be checked by 
computing the pseudo probabilities for different components, eq. (1311) . at several values of 
the binding energy. To this end, we apply the transformation eq. (1251) (or eq. fT26|) for the 
3 Si — 3 D\ state) to the set of functions g n , preliminarily obtained numerically from the 
series (l35l) - (l3~Tj) (or (|38|) - (|4T1) ). The corresponding results for the 1 So channel are collected 
in Table IXlTl where the pseudo probabilities for different components have been computed 
assuming the total normalization to be 

P++ + P— +P e + P Q = 1. (63) 

As expected, the contribution of the 1 Sq~ component is negligibly small for weakly bound 
systems and significant for large binding energies. Similar conclusions hold for the 3 S\- 
3 Di state as well. As mentioned, the BS system of equations for the partial components 
have been solved for a coupling constant g below the critical value. However, there is also 
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a lower limit below which the solution of the BS equation disappears, in a fully analogy 
with the nonrelativistic Schrddinger equation, when, for shallow potentials, bound states 
do not exist. Another observation is the result that the condition det(g 2 A — 1) = 
provides an equation for M(g 2 ) which in the interval < M < 2m can have more than 
one solution. In such a case the lowest value of M corresponds to the ground state, while 
others refer to the discrete excited state of the system. We investigated also the behaviour 
of the partial vertex functions g\, . . . , for the excited states and found that, likewise 
in the nonrelativistic case, these functions possess zeros as functions of This is in a 



good agreement with our previous results obtained for the two-dimensional mesh [33]. 



D. Pseudoscalar meson exchange 



As follows from eq. (1321) . the system of equations for pseudoscalar meson exchange is 
quite similar to the previously described scalar case (see also Tables HVTlVIl) . However, 
contrary to the scalar case, when the main components g\ and #2 in eqs. (1421) and (|4"5T) 
are determined by pure attractive kernels, for the pseudoscalar exchange one of these 
components is governed by a pure repulsive kernel (c.f. Tables HT1 and [TTT]) . Thus, the 
resulting balance of forces forming the bound state for pseudoscalar exchange is more 
sophisticated. As a consequence, to ensure an attractive residual kernel for creation of 
a given bound state M, the BS equation requires larger values of the coupling constant 
g 2 . Consequently, this can lead to subtle situations when the minimal value of g 2 is close 
to or even above the critical value g 2 r , which may cause problems with the stability of 
the solution. Namely such a situation we encountered in our numerical calculations for 
a two-nucleon like system with pseudoscalar exchanges for which the stable bound state 
(without cutoff form factors) does not exist at all. Similar result has been obtained also 
within the LF approach [8| and seems to be of a general nature, i.e. there is no relativistic 
bound state in a deuteron-like system with pure pseudoscalar exchange. 

Another striking feature here is connected to the behavior of the binding energy as a 
function of the coupling constant. In Fig. [HI the dependencies of the binding energy for 
the 1 5'o state are shown for some values of the cutoff parameter A. It is clear that such 
a dependence is very sharp - binding energyjncreases very rapidly with g 2 , which is in 
agreement with the results reported in Ref. 
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E. Vector meson exchange 



Since the contribution of vector meson exchange in the nucleon-nucleon potential is 
repulsive a deuteron-like bound state cannot be formed by pure vector meson exchanges, 
therefore an investigation of the homogeneous BS equation with such kernels is hampered. 
However, one can consider a different two-fermion system such as the electron-positron 
pair, for which the vector exchange potential does have a bound state. Note that the 



general form of the BS equation, eq. (TTTT) . in t 



almost unchanged (see, e.g. Refs. 



re particle-antiparticle channel is maintained 



3211 and 



41] ) so that a relativistic description of the 



positronium can be achieved by eq. (ITT]) with T(l) = 7^, f(2) = •y u and V(p, k) = 
9nv/{p ~~ k) 2 - Unfortunately, the knowledge of the positronium, as a two-fermion bound 
system, is rather scarce [20| , even the absolute value of the binding energy is not definitely 
established. 

Nowadays only the transition energy between different positronium levels (e.g., be- 
tween para- and ortho-positronium) is an object of experimental and theoretical inves- 



tigation 



20J. Compared to the electron mass this quantity is very small 21|, being of 



the order AB ~ m e a® m /A ~ lCT 13 m e (where a em is the fine structure constant, m e is 
the electron mass and B ~ m e «g m /4 is the positronium binding energy predicted by 
the nonrelativistic Schrodinger equation). Consequently, the procedure of solving the BS 
equation numerically for the bound state M ~ 2m e — B and calculations of the transition 
energy AB and comparison with experimental data require extremely high accuracy of 
calculations and large computational resources. Moreover, as mentioned before, in case 
of vanishing masses of the exchange particle, the convergence of the method is rather ill 
defined and the implementation of cut-off form factors becomes a necessity in numerical 
calculations. Therefore, in analyzing effects of relativistic corrections computed within 
different schemes (e.g., within perturbative quantum electro-dynamics (QED), LF dynam- 
ics and BS formalism) one usually solves the corresponding equation with cut-off form 
factors and at an effective coupling constant g much larger than the fine structure constant 
a em (see Ref. js|) and compares the obtained results with those known analytically from 
the nonrelativistic Schrodinger equation. For this reason, a commonly accepted value 
of g is g 2 = 3.77, which corresponds to a — j- — 0.3 and to a nonrelativistic binding 
energy B = m e a 2 /A = 2.25 • 10~ 2 m e [8|. In our calculations we also adopted such a 
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coupling constant g 2 = 3.77 for which the para- and ortho- positronium binding energy 
has been calculated. The results are presented in Table IXHIt where the cut-off param- 
eter A and binding energies are given in units of the electron mass m e . The obtained 
results demonstrate that at large enough values of A, the solution of the BS equation is 
quite stable and almost independent on A. It is also seen that the BS equation provides 
a binding energy almost twice smaller than the nonrelativistic one which implies that 
the relativistic corrections are of repulsive nature. Analogous result has been obtained 
within the LF dynamics for the 0~ positronium j^ . However, calculations within the 
perturbative QED show that the first order relativistic corrections are attractive. This 
contradictory result can serve as an indication that the adopted interaction kernel is not 
accurate enough for a refined description of the positronium within the ladder approxima- 
tion. Other channels (e.g, the electron-positron annihilation 
the ladder approximation should be included into the analysis. 



411 ]) and/or terms beyond 



V. CONCLUSION 



We generalize a method based on hyperspherical harmonics to solve the homogeneous 
spinor-spinor Bethe-Salpeter equation in Euclidean space. To do so, we introduce a new 
basis of spin-angular harmonics, suitable to expand the Bethe-Salpeter vertex into four- 
dimensional hyperspherical harmonics. We obtain an explicit form of the corresponding 
system of one-dimensional integral equations for the partial components and formulate a 
proper numerical algorithm to solve this system of equations. The BS vertex functions 
are studied in detail for the 1 So and 3 Si — 3 Di bound states with scalar, pseudoscalar and 
vector meson exchanges. Our results are in a good agreement with calculations within 
the non relativistic and Light Front Dynamics approaches. 

Within the novel method the effectiveness of the numerical procedure is analyzed for 
the scalar, pseudoscalar and vector meson exchanges and conditions for stability of the 
solution are established. It is demonstrated that above some critical values of the coupling 
constant the solution of the BS equation does not exist unless the cut-off form factors are 
considered. 

An advantage of the method is the possibility to present the numerical solution in 
a reliable and simple analytical parameterized form, extremely convenient in practical 
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calculations of matrix elements within the BS formalism and for analytical continuation 
of the solution back to Minkowski space. 

The method allows us to describe, in a covariant way, realistic two-body systems, such 
as the deuteron, the positronium and the variety of known mesons, as bound states of 
quark- ant iquark pairs and to solve the inhomogeneous BS equation for scattering states 
in the continuum. 
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APPENDIX A: PARTIAL KERNELS 



Here we present the explicit form of the kernels V k , k (p) which determine the partial 
kernels A^(p) and B^{p), eqs. fj42l - fH4l) and fH5l) - fTTfj) . For this let us introduce an 
auxiliary quantity S l k , k (p) defined as 



7T TT 



sin 2 xX k >i{x) Xki{x) 

\2 







(p 2 + m 2 - M 2 /4y + M 2 p 2 cos 2 x 



2 /! (-2) m / (£;' + l)(fc'-Q! (fc + l)(fc-Q! 

npM(p 2 + m 2 - M 2 /4)y (k' + I + 1)! (jfe + Z + l)! 

x C l +Uiz)(z* + l)^Q^ +| (zz), (Al) 

where max (min) is the maximum (minimum) index of fc, fc', C^^^iz) and Q^^ii^) 
stand for the Gegenbauer polynomials and Legendre functions of the second kind, respec- 
tively, of imaginary argument iz with 

_ p 2 + m 2 - M 2 /4 
Z ~ pM ' 



28 



Note that, as follows from properties of Gegenbauer polynomials, the quantities S l k , k are 
different from zero for nonnegative k, k' and k + k' integer. The Legendre functions in eq. 
(I AID for / = 0, 1, 2 explicitly read as 



i 

I 2 

'"+3 
3 



Ql +3 _(iz) = (n + l)zQl^_(iz) - z(n + 2)Ql^_(iz) 



1 



Ql + ,{iz) = (n + l)zQl^(iz) - i(n + i)Ql^(iz) 



' n+ 



l 
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By using the recurrent relations for the Gegenbauer polynomials 39|, the partial kernels 
V k , k can be expressed via S l k , k as 



^13 - /WvMM^v- 1 J+ 2 o0 1 / 3 q0 



l y/(j-l)(j+3) / /j + 2 / j \ 1 / {j - l)j c0 

4 + 3 + Vj + 2 y l^ + 4\/(j + 2)( J + 3)^ +2 ' 

t/22 = f d ■ a ^fco(x)^o(x) 2 _ 1 / cO , oO , oO , 9 o \ 

^fcj — y ax sin X^ 2 + ^2 cog 2 y cos ^ - + + D fc+i,i-i + D fe+ij+iJ ; 

V kj = / 4sm X -72~i P2 r~ cos X sm X 

J J A 1 + B z cos z x 

_ If / + ^ r.0 _ / cO j_ . / + ^ qO _ / r-0 1 

^ - / ^sin 2 X ^f^ x oos X = 1 - ( Cl (k)Sl +ld + c 2 {k)Sl_ Xj ) , 

iff S /^sin 2 x ^f j2( 2 X) sin X 

J J A z + B z cos z x 

_ 1 / (j + 2)(j + 3) cl _ 1 / (j - l)j ! 

" 2 y + 2V(j + l)(j + 2)^ +1 ' 

^ - / rfxsin 2 x j^^ cos 2 x = J Hk) Cl (j)Sl +l!j+1 + c.W^O-)^-! 

Vf - / ^»in 2 X ^ M + ( ^f^ oos 2 x= ^ (di(fc)diC7)g fc Wi + di(fc)d 2 (j-)gWi 

+ c?2(A;)c?i(i)^fe_ij+i + d 2 (k)d 2 (j)Sl_ hj _ 1 ) , 

t,57 /" . -2 ^fcl(x)^'2(x) 

^fc, = / «xsm x-p^ n9 9 cos x sin x 

3 J A 2 + B 2 cos 2 x 



1 m / (j + 2)(j + 3) x / 



1 m / (j + 2)( J + 3) (j-l)j 

i C2W V i(JTi) Vl - 1 "V0' + i)0- + 2)' 



, fc(fc + 3) / (fc-l)(fc + 2) 

Cl(A ° = ^F#T2)' C2(A;) = V k(k + i) ■ 
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M 2 /4, 
■ l)(fc + 4) 



(A: + l)(fc + 2)' 



B = pM, 
d 2 (k) 



(fc-2)(fc + 3) 
+ 1) ' 
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'So 


1 C++ 1 c 3 pe 3 po 

Dq Dq fg -Tq 


3 S 1 - 3 D 1 


3 5++ 3 sf- 3 L>+ + 3 L>" 

3 pe 3 po 1 pe 1 po 



TABLE I: Classification of the p-spin partial components for the 1 So and 3 Si~ 3 Di channels in 
spectroscopic notation. 



n 


b n [S] 


b n [PS] 


bn[V] 


1 


1 


-1 


4 


2 


1 


1 


-2 


3 


1 


-1 





4 


1 


1 


-2 



TABLE II: The coefficients b n , eq. (f32j) . in the 1 So channel with scalar (S), pseudoscalar (PS) 
and vector (V, without the tensor part) meson exchanges. 
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n 


j r m 

K[S] 


b n [PS] 


6n|V] 


1 


l 


1 


2 


2 


l 


-1 





3 


l 


1 


-4 


4 


l 


-1 


-2 


5 


l 


1 


2 


6 


l 


-1 





7 


l 


1 


2 


8 


l 


-1 






TABLE III: The same as in Table U but for the 3 S 1 - 3 D 1 channel. 



n 


a n i 


a n 2 


«n3 


Un4 


1 


Ml + m 2 + ~k 2 


mM 


-M\k\ 





2 


mM 




-2m\k\ 


-2k A \k\ 


3 


M\k\ 


2m\k\ 




2mk A 


4 





-2k 4 \k\ 


—2mki 





TABLE IV: The quantities a nm , eq. for the 1 5'o state. 
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n 




«n2 


On3 




1 




mM 




M\k\ 


2 


mM 


m2 _L TO 2 i 4r,2 A: 2 
— -f m -f g«, 4 — 


M|fe| 


-^fm\k\ 


3 




^M|fc| 







4 


fM\k\ 


2^m|fc| 







5 







2^m|fc| 


-^M|fc| 


6 





-^|fc| 2 


^M\k\ 


-^fm\k\ 


7 












8 










2m/c4 



TABLE V: The quantities a nm , eq. (f34l) . for the 3 Si~ 3 Di state. 
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n 


«n5 


«n6 




On8 


1 












2 












3 


-2^m|fc| 


M|fe| 


— 2mfc 4 





4 


f M|fc| 







—2771^4 


5 


^f- + m 2 + |/c| — 









6 


mM 


Ml+m 2 + + f 







7 


-^h\k\ 










8 










_ + m 2 _ 2 Jfe2 + jfc2 



TABLE VI: Continuation of Table El 
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n 


Anl 


4n2 


An3 
A jj' 


AnA 


1 


-{M 2 /A + m 2 + k 2 ) 


-mM 


kM 







X ^2j-2,2j'-2 


x 5 , 2j-2,2j'-2 


* V 2j-2,2j'-l 








-(M 2 /4 + m 2 -P) 






2 




x ^2j-2,2j'~2 
v 2j-2,2j'-2 


2mk V2j-2,2j'-i 


ZK V 2j',2j-2 


3 


-kM 


—2mk 


(M 2 /4 - m 2 + fc 2 ) 


—2mk 






x v 2j'-2,2j-l 


XD 2j-l,2j'-l 


xV 34 

X V 2j-l,2j> 


4 





Ik 2 V 24 

ZK V 2j,2j'-2 


2/A; V$ j2jl _ x 


(M 2 /4 - m 2 - k 2 ) 

X °2j,2j' + ZK K 2j',2m 



TABLE VII: The partial kernels A^'(fc) denned by eqs. (02} -(@U). For the explicit form of the 
introduced quantities S^, k (k) and VJj?' (k) see Appendix lAl 
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n 


33 


J/ 


5"^ 


B nA , 

33 




-(M 2 /4 + m 2 + P/3) 




2mA; 


V2Mk 


1 


A °2jf-2,2j'-2 

— 2&; 2 V 2 22 _ 2j2? -,_ 2 /3 


-mMS^_ 22j ,_ 2 


x^ 3 - 2)2 ,-'-i/V3 


x^ 3 -W-i/V3 






-(M 2 /4 + m 2 - fc 2 /3) 


MA; 


2\/6 m/c 


2 




X ^2j-2,2j'-2 

-4fc 2 Ff_ w _ 2 /3 


x^- 2A -'-i/>/3 


X ^2 1 /-2,2j'-l/ 3 


3 


— 2mfc 




(M 2 /4 - m 2 + P) 









X^-2,2,-l/^3 


X^j-l^j'-l 




4 









(M 2 /4 - m 2 + fc 2 ) 




X ^2/-2,2j-l/ 3 


xV 2 1 ^_ 22j _ 1 /3 




X ^2j-l,2i'-l 


5 


-2v / 2fc 2 y 2 1 /„ 2>2j /3 





— 2\/6 m& 

X ^2^-l,2j/ 3 


fcM^_ lj2i /V3 


6 





2V2*V$_ a , 2j /3 


X ^2^-l,2i/ 3 


2m^f_ 12j /^3 


7 


-2py 2 ^.,_ 2 /V3 





-2mfcy 2 ^,_ 1 





8 





-2V6fc 2 ^4 2j; _ 2 /3 





-2mkV^ 2jl _ x 



TABLE VIII: The partial kernels Bp ' (k) denned by eqs. (@5])-(|I7|). For the explicit form of the 
quantities S^ : -,(k) and VJ*?' (k) see Appendix lAl 
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n 


5 n! ! 

33 


33 




jr 


1 


-2V2~k 2 V 2 \ 5 _ 2<2m /3 





-2p 

X^2?,2i-2/^3 





2 





2^2py 2 iL 2 , 2i /3 





-2-Sk 2 

X ^2j,2j'-2/3 


3 


2V6mA:y 2 3 /_ li2j ,/3 


V6fcMy 2 3 /_ lj2j ,/3 







4 


-fcMy 2 3 /_ li2i ,/V3 


-2mkV£f_ h2jl /V3 







5 


-(M 2 /4 + m 2 -P/3) 
X^. 2i ,-4p^ 2 .,/3 


-mMS\ h2y 


x^.2,/3 





6 




-(M 2 /4 + m 2 + jfc 2 /3) 

X ^2j,2j' ~ 2^ 2 ^2j^2j'/3 





-2fc 2 


7 


2x/6fc 2 F 2 ^.,/3 





(M 2 /4-m 2 -fc 2 ) 

XS 2j,2j' 

+2k 2 V 2j ^ 2j , 





8 










(M 2 /4 - m 2 - fc 2 ) 

X5 2j,2j' 



TABLE IX: Continuation of Table |Vffll 
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9 -i r 

g 2 = 15 


/x = 0.15 GeF/c 2 


jj, = 0.5 GeV/c 2 




iV G = 32 


N G = 64 


N G = 96 


N G = 32 


N G = 64 


N G = 96 


1 


1.9399 


1.9399 


1.9399 


1.9984 


1.9984 


1.9984 


2 


1.9370 


1.9370 


1.9370 


1.9982 


1.9982 


1.9982 


3 


1.9368 


1.9368 


1.9368 


1.9982 


1.9982 


1.9982 


4 


1.9368 


1.9368 


1.9368 


1.9982 


1.9982 


1.9982 




g 2 = 30 


/x = 0.15 GeV/c 2 


fi = 0.5 GeV/c 2 


Mnax 


N G = 32 


N G = 64 


N G = 96 


N G = 32 


N G = 64 


N G = 96 


1 


1.7932 


1.7910 


1.7905 


1.9167 


1.9142 


1.9137 


2 


1.7897 


1.7875 


1.7871 


1.9152 


1.9127 


1.9122 


3 


1.7896 


1.7874 


1.7870 


1.9152 


1.9127 


1.9122 


4 


1.7896 


1.7874 


1.7870 


1.9152 


1.9127 


1.9122 



TABLE X: Convergence of the calculated values of mass M with respect to the number of points 
in Gaussian mesh N G and the number of hyperspherical components M max for the values of the 
coupling constant g 2 = 15 and g 2 = 30 in 1 Sq channel for the scalar meson exchange kernel. 



j 




4 




4 


a\ 


x 2 


1 


0.8162 


3.7476 


-3.3910 


2.9564 


-0.0107 


3.4 


• 10- 4 


2 


0.5875 


-3.9733 


6.6483 


-6.5271 


3.6482 


1.5 


• 10~ 6 


3 


0.4556 


4.2805 


-5.6270 


6.1247 


-4.8948 


2.4 


• 10~ 6 


4 


0.3680 


-4.2026 


9.4868 


-22.191 


17.300 


8 • 


10~ 5 



TABLE XI: Numerical values of the parameters in formula (|6ip and the corresponding x 2 values. 
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9 2 


M, GeV/c 2 


P++ 








Po 


Pe 


15 


1.937 


1.012 


-1.18 


10~ 


-3 


-6.63 • 


io- 


-3 


-4.37- 


io- 


-3 


20 


1.892 


1.020 


-2.99 


10~ 


-3 


-1.07- 


io- 


-2 


-6.92 • 


io- 


-3 


25 


1.842 


1.030 


-6.22 


io- 


-3 


-1.46 • 


io- 


-2 


-9.41 • 


io- 


-3 


27 


1.820 


1.034 


-8.11 


io- 


-3 


-1.61 • 


io- 


-2 


-1.03- 


io- 


-2 


29 


1.798 


1.039 


-1.05 


io- 


-2 


-1.75- 


io- 


-2 


-1.12 • 


io- 


-2 


29.5 


1.788 


1.041 


-1.25 


10- 


-2 


-1.80 • 


io- 


-2 


-1.16 • 


io- 


-2 


29.6 


1.5 


1.210 


-0.19 




-1.24- 


10- 


-2 


-7.77- 


10- 


-.3 



TABLE XII: Pseudo-probabilities, eq. (|3ip . in the 1 Sq state at given g 2 and M, normalized 
as (|63p . In the LF formalism only " + +" and "o" components are relevant to the computed 
pseudo-probabilities . 



A, m 


c 2 


10 


20 


30 


40 


50 


60 


70 


80 


Bl s ■ 


10 2 


1.378 


1.362 


1.360 


1.355 


1.352 


1.350 


1.349 


1.348 




10 2 


1.273 


1.268 


1.262 


1.260 


1.259 


1.258 


1.257 


1.256 



TABLE XIII: Binding energies of positronium states as a function of a cutoff parameter A in 
units of m e c 2 at g 2 = 3.77. 
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0.0 0.5 1.0 1.5 2.0~ 2.5 3.0 
p [GeV/c] 



FIG. 1: Functions g\,j = 1,...,4, eq. ([35]) . at the gaussian mesh points Nq = 96. Closed 
squares correspond to g\, closed circles - g\ multiplied by 10, triangles -g\ multiplied by 100, 
open circles -g\ multiplied by 1000, the solid lines correspond to the fitted functions g\ by 
formula (|6ip . The overall normalization constant is arbitrary. 
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FIG. 2: Mass of the bound state in the 1 S'o channel as a functions of the cutoff parameter A for 
different values of the coupling constant g at m = 1 GeV/c 2 , n = 0.15 GeV/c 2 
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1 2 3 4 5 
p [GeV/c] 

FIG. 3: Hyperspherical components g{,j = 1,2,3, for A = 5 GeV/c and coupling constant 
g 2 = 30 (solid line) and g 2 = 48 (dashed line). For the component g\ the lines reproduce the 
results of interpolation by eq. ([62]) . 
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FIG. 4: Masses of bound state in the 3 5i- 3 .Di channel as functions of the cutoff parameter A 
for different values of the coupling constant g. Critical value g^ T ~ 78. 
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FIG. 5: Dependence of the mass of the bound state M in the 1 Sq channel for different approaches. 
The solid line represents results within the BS formalism, dashed line corresponds to Light 
Front (LF) dynamics and the dotted line is the result of the non relativistic (NR) Schroedinger 
approach. 
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FIG. 6: Comparison of the BS partial "++"-component eq. (|27p. solid line, with LF partial 



component /i Ref. 



4C], dashed line, at the binding energy B = 2m — M = 1 MeV/c 2 
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FIG. 8: Sensitivity of the binding energy B = 2m — M to the coupling constant g 2 and to the 
cutoff parameter A for the 1 Sq state with pseudoscalar meson exchange. 
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